Fronts of spin tunneling in molecular magnets 
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Dissipative spin-tunneling transitions at biased resonances in molecular magnets such as Mni2 
Ac are controlled by the dipolar field that can bring the system on and off resonance. It is shown 
that this leads to spin relaxation in form of propagating fronts of tunneling, with the dipolar field 
adjusting self-consistently to provide a zero bias within the front core. There are two regimes of 
the front propagation: laminar and non-laminar with discontinuous magnetization and dipolar field. 
In the laminar regime the speed of the front can exceed that of the magnetic deflagration, if the 
transverse field is large enough. Fronts of tunneling can be initiated by magnetic field sweep near 
the end of the crystal. 

PACS numbers: 75.50.Xx,75.45.+j,76.20.-fq 
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I. INTRODUCTION 

Molecular magnets (MM), including their first and 
mostly studued representative Mni2 Acj^ have initially 
attracted attention as molecules vsfith the effective big 
spin 5* = 10 showing bistability as a result of a strong uni- 
axial anisotropy —DS^.^ Resonance spin tunneling man- 
ifested in the magnetic hysteresis loops''-''-^ with steps at 
the field values B Bk = kD/ (g^ig) , fc = 0, ±1, ±2, . . . 
made molecular magnets a hotspot of research during 
more than 10 years. 

Crystals of molecular magnets do not show a significant 
exchange interaction because the magnetic core of the 
molecule is screened by organic ligands. Thus magnetic 
molecules remain largely superparamagnetic, although 
MM can order below 1 K due to dipole-dipole interac- 
tions (DDI),^ii 

An important role of the DDI is that the dipolar field 
created by the spins is large enough to change the res- 
onance condition for the up and down spins and thus 
to strongly influence spin tunneling. Fully ordered spins 
in an elongated Mni2 Ac crystal create the dipolar field 
~ 52.6 mT at a molecule<ii^ This becomes com- 
parable with the resonance width defined by the tunnel 
splitting A in transverse magnetic fields above 5 T, for 
the k = 1 tunneling resonance. For smaller transverse 
fields, A is much smaller and thus the DDI can com- 
pletely block the resonant tunneling. The action of the 
dipolar field is dynamical and self-consistent since tun- 
neling of spins causes the dipolar field to change, blocking 
or allowing resonant transitions. 

The role of the DDI in spin tunneling was recognized in 
Refs. 9.10,11,12,13,14,15,16 where Monte Carlo simula- 
tions were done on the basis of a phenomenological model 
involving discrete jumps of the spins through the instan- 
taneous "tunneling window" . The main purpose of these 
studies was to explain the \/t relaxation experimentally 
observed in Mni2 Ac.^^ 

It was not understood until recently that DDI in 
molecular magnets can result in spatially-inhomogeneous 
states creating the dipolar field such that the system is on 
resonance in some regions of space where spins can relax, 
leading to moving fronts. An example is the domain wall 



in elongated dipolar-ordered crystals of Mni2 Ac at low 
temperatures. The reduced dipolar field at T = in Fig. 
2 of Ref. .7 is close to zero in the region around the cen- 
ter of the domain wall with the width of the order of the 
crystal's thickness. It should be mentioned that Mni2 Ac 
remains the only molecular magnet that can be grown in 
long crystals required for such kind of phenomena. 

Similar effects can take place in the tunneling at bi- 
ased resonances, A; > 1. If the external field approaches 
the resonance by a slow sweep, as was the case in many 
experiments, moving walls of tunneling can be created 
near the ends of long crystals (where the dipolar bias is 
smaller) and then penetrate into their depth with a speed 
unrelated to the sweep ratej^^ The role of the sweep is 
only to create an initial state for the wall of tunneling 
to start. It was argued that this mechanism can explain 
the width of the steps in dynamic hysteresis curves^*^ 
by the time needed for the wall of tunneling to cross the 
crystal. Non-uniformity of the magnetization in Mni2 
Ac developing during spin tunneling was detected by lo- 
cal measurements earlier.— 

The walls or self-organized patterns of spin tunnel- 
ing investigated in Ref. [H are not exactly fronts be- 
cause they are lacking the combined space-time depen- 
dence on the argument z vt only, where v is the speed 
of the front. Frozen-in quasiperiodic spatial structures 
have been found behind these moving walls. In fact, true 
smooth fronts of spin tunneling do exist in the range of 
the external bias smaller than that in Ref. Ts'. Study- 
ing these fronts and their transition to the moving walls 
with a nonuniformity behind with increasing the bias is 
the purpose of this article. It will be shown that there 
are two regimes. 

For the external bias not exceeding a critical value, the 
true fronts (that can be called "laminar" ) are realized in 
which the dipolar field adjusts to create a resonance in 
the front region. In the limit of strong dipolar field (rel- 
ative to the resonance width) the front speed and the 
magnetization behind the front can be calculated analyt- 
ically and are independent of the strength of the DDI. 

For a larger external bias, the magnetization distribu- 
tion and thus the dipolar field in the wall cannot fully 
adjust to provide the resonance condition. In this case 
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the wall is moving with a quasiperiodically varying speed 
leaving a quasiperiodic state behind. The average wall 
speed decreases with the DDI strength quadratically. 

The dipolar mechanism of spin tunneling is resembling 
magnetic deflagration in Mni2 Ac.^"-^^ Here, instead of 
the temperature, the relaxation rate is controlled by the 
self-consistent dipolar field bringing the system on or off 
resonance. Thus, in a sense, one can call the phenomenon 
studied here cold deflagration. Of course, the heat release 
in the course of the cold deflagration can give rise to the 
regular deflagration, especially for high resonances k and 
well thermally isolated crystals. In this case the two kinds 
of deflagration can compete. 

The rest of the article is organized as follows. In Sec. 
|ll]the dynamics of spin tunneling between the metastable 
ground state and a resonant excited state on the other 
side of the barrier is considered. The simplified over- 
damped equations of motion are obtained in the case 
of the tunnel splitting frequency A/h smaller than the 
damping of the excited state F. It is further argued that 
in the presence of disorder that spreads resonances one 
can use overdamped equations in a generalized form also 
for larger A. In Sec. Illll the dipolar field created by a wall 
of magnetization is calculated for the cylindrical and rib- 
bon geometries. In Sec. |IV]the full system of cold defla- 
gration equations is written and transformed into dinien- 
sionless form. In Sec. |V]the limit of strong DDI is studied 
and analytical expressions for the residual magnetization 
behind the front and the front speed are obtained. Sec. 
I VII provides the results of numerical calculations in both 
regimes of the wall propagation. 

II. SPIN TUNNELING AND RELAXATION 

We will be using the generic giant-spin model of molec- 
ular magnets with the Hamiltonian 

H = -DS^ - g^igB^Sz - g^i^B^Sx + ■ ■ ■ , (f) 

where D is the uniaxial anosotropy and 

B - Be,t + B(^) (2) 

is the total magnetic field, including the external and 
dipolar fields. Suppressed terms in the Hamiltonian can 
include the biaxial and fourth-order anisotropy that can 
make a contribution into the tunnel splitting A of the 
resonant spin up- and down states. Since the most inter- 
esting situation arises in the case of a large A that can 
only be created by a strong transverse field, the dropped 
terms will not be needed. For B^: = the exact quantum 
states of H are |m) with —S<m<S, their energies 
being = —Dm? — g^gBzin. The resonance condition 
£m = £m' between all states |m) on the left side of the 
barrier (m < 0) and |m') on the right side of the barrier 
{m' = — m — k) is satisfied for the resonance fields 

B.^Bk, Bk = kD/{gfiB), A: = 0,1,... (3) 

This resonance condition turns out to be independent of 
the transverse field. 



Application of the transverse field leads to the two ef- 
fects. First, each state |m) hybridyzes with neighbouring 
states within the same well forming the state that can 
be denoted as IV'm) • Physically this corresponds to spin 
canting in the direction of the transverse field. Second, 
the states lip^) on different sides of the barrier hybridize 
because of the resonance spin tunneling near B^ = Bk ■ Of 
course, one can speak of the states \ip^) for not too strong 
transverse field, so that there still are low-lying states well 
localized within one of the wells. The states IV'm) P'"0- 
vide a basis for a simplified treatment of spin tunneling 
and relaxation near resonances that otherwise has to be 
considered within the density-matrix formalismj^i^ 

Consider the metastable ground state of a molecular 
magnet, |V'm) = |V'-s) ' near a tunneling resonance with 
an excited state [ipm/) on the right side of the barrier. 
Fig. I of Ref. IS. The dynamics of tunneling at low 
temperatures is described by the subset of the density 
matrix equation (DME) taking into account only these 
two levels. The level IV'm') '^^^ decay into lower-lying 
levels within the same well with rate Tm'- Since there 
are no incoming relaxation processes for the state \ipm') 
at low temperatures, the DME can be simplified to the 
form of the damped Schrodinger equation 

i A 

C-S - -2^Cra' 

c^' = [—-2^^n')c„,---c^s- (4) 

Here A is the tunnel splitting and W = e„i — Em' is the 
energy bias between the two levels, 

W={S + m') g^ls (b, - B,, + B^^^) = W,^, + W(^). 

(5) 

In fact, here one should use the values of S and m' cor- 
rected for spin canting. Since j^/i-s) is the lowest state 
in the left well, it cannot decay. The numbers of particles 
in the states are defined by 

n^s = \c-s\^ , rim' = |c„i'|^ (6) 

etc. The spin polarization in our low-temperature tun- 

nehng process is given by {S^.) = -Sn_s + Y.r,i=m' mn,n- 
As the states with m — m' + 1, . . . , S~l decay faster than 
\4'm') ' their contribution can be neglected. Then for the 
normalized spin-average variable 

= {S.) /S, (7) 

one obtains 

a, = l- 2n_s - (1 - m/S)nra' ■ (8) 

In the overdamped case Tm' 3> A/ft the variable c„i' in 
Eq. ([4]) adiabatically adjusts to the instantaneous value 
of c_s. Setting Cm' = in the second of these equations 
one obtains 



2hWlh + iTm'/'2.' 
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Inserting it into the first of equations ^ one obtains 
iA^ c-s 



c-s 



Ah? W/h + iVm'/'^' 



(10) 



With the help of Eq. one obtains the equation for 
the metastable population n^s in the form 



n = — Fn, 



(11) 



where the subscript ~S has been dropped for trans- 
parency and the dissipative tunneling rate V is given by22 



A2 



(12) 



In the overdamped limit one has Um' *C n-s = n, so that 
Eq. ([5]) simplifies to 



(7, = 1 — 2n. 



(13) 



The decay rate Tm' is mainly due to the relaxation 
between the adjacent energy levels in the right well, 
Tm' = r^'.m'+i, see Eq. (A9) of Ref. H or Eq. (294) 
of Ref. f23|. For fc = 1 and thus m' = 9 one has T^, w lO"^ 
s~^, while A/h reaches a comparable value in the trans- 
verse field above 3 T. At higher transverse fields the tun- 
neling dynamics should be underdamped. Resonances 
with higher k have larger splitting A and become under- 
damped in smaller transverse fields. 

In the case of underdamped resonances, A/ft > r„i', 
the rate of dissipative spin tunneling can be described by 
the integral relaxation time Tint resulting in the effective 
rati 
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A2 



Tint 2ft2 ^2 + {T^,/2f 



where 



A^ 



(14) 



(15) 



One can see that in the underdamped case the width of 
the Lorentzian becomes of the order of (A/2) /ft, com- 
pared to Vm' /2 in the overdamped case. Although Eq. 
PT|) with r given by Eq. ([T3|) does not accurately de- 
scribe the oscillating dynamics of the system in the un- 
derdamped case, in particular the Landau-Zener effect, 
it will be used below as an approximation for the many- 
body problem with coupling via the dipolar field in both 
overdamped and underdamped cases. A more rigorous 
approach based on Eq. ^ requires much more computer 
time because of fast oscillations. On the other hand, 
oscillations at tunneling resonances have never been ex- 
perimentally observed in MM because of the resonance 
spread as a result of ligand disorder and other factors. 
For the low-bias resonances such as fc = 1 and thus 
77i' = — 1 the contribution of 7i„i' in Eq. ([8]) can be 
neglected, thus Eq. (fT5|) will be used in all cases. 



III. DIPOLAR FIELD 

The dipolar field and ensuing dipolar bias of tunneling 
resonances in crystals of molecular magnets have been 
discussed in detailes in Ref. 0, so that only a short sum- 
mary with necessary changes will be provided below. 

The z component of dipolar field at site i (i.e., at a par- 
ticular magnetic molecule), created by molecular spins 
polarized along the z axis is given by 



B. 



vo 



■D, 



where vq is the unit-cell volume and 



= vo- 



3(e^-ny) -1 



(16) 



Z^- (17) 



Inside a uniformly magnetized ellipsoid the dipolar field 



is uniform and one has 

5„ = 51?'^) + 



Dzz<7z, where 



Am/ 



(l/3 



n 



(2) 



(18) 



v is the number of magnetic molecules per unit cell and 
n(^) = 0, 1/3, and 1 for a cyhnder, sphere, and disc, 

respectively. D'fz'^^ depends on the lattice structure. For 
Mni2 Ac lattice summation yields D^zz^'' = 2.155 that 
results in Z)i7'^ = 10.53 for a cylinder. Then Eq. HID 

yields the dipolar field ^ 52.6 mT in an elongated 
sample that was also obtained experimentally.^ 

The dipolar bias VF^^^ in Eq. can be written in the 
form 



D) 



1 



m 
'J 



EdDzz, 



where 



Er 



iSg^J■B) ho 



(19) 



(20) 



is the dipolar energy, Eo/kB = 0.0671 K for Mni2 Ac. 
Since the dipolar field depends on the magnetization and 
its values in an elongated crystal can change between 
—52.6 mT and 52.6 mT, one can conclude that, accord- 
ing to Eq. (O, the resonance condition W — can be, 
in principle, satisfied in the dipolar tunneling window 
around the resonance —52.6 mT < Bz — Bzk < 52.6 mT. 
This dipolar window is much smaller than the distance 
between the two tunneling resonances that is about 0.5 
T. Practically, for a negative external bias Bz — Bzk the 
relaxation is hindered by the causality: To produce a 
positive dipolar field that would balance the negative ex- 
ternal bias, spins should already be on the right side of 
the barrier. 

For a cylinder of length L and radius R with the sym- 
metry axis z along the easy axis, magnetized with ct^ = 
(7z{z), the reduced z-field along the symmetry axis has 
the form 

f^^^ , / 2'KvR^(Jz{z') , , , , 

Dzz{z)= I dz'- ^ -kaz{z), (21) 



-L/2 



(z'-z)'+i?2 
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where 



k EE 87ri//3 - D'^^p^^ = Attv - i:*^^') > 0, (22) 

fc = 14.6 for Mni2 Ac. In particular, for a long uniformly 
polarized cylinder one obtains 

a, + (^i^*^) - 2 (7,, (23) 



2-KVZ 
+ i?2 



At the left end of the long cylinder, z = 0, one has 
Dzz — [d'^zF^^ ~ 2'Kv/'i^ Oz. For Mni2 Ac one obtains 



Dzz = — 2.03cr2, having the sign opposite to that of the 
field in the depth, D^zl^'' Oz- This means that a homoge- 
neously magnetized state in the resonance external field 
Bz = Bzk may be unstable with respect to spin tunnel- 
ing beginning in the vicinity of the ends of the cylinder, 
as at some point near the end the resonance condition 
Dzz{z) = is satisfied. To the contrary, dipolar field in 
the depth of a uniformly magnetized cylinder provides 
the dipolar bias that puts the system off resonance and 
makes the transition rate T very small. If the external 
field Bz is swept in the positive direction towards the 
resonance, spin tunneling begins near the ends of the 
crystal and then it propagates inside the crystal as a 
moving wall of tunneling, as the dipolar field changes 
self-consistently.^® 

One can also calculate the dipolar field on the z sym- 
metry axis of a slab of length L and sides a and b. In the 
case a ~ 6 the results are similar to those for the cylinder 
above. For a thick slab (a ribbon) with a <C 6 one obtains 



Dzz{z) = / dz' — 

J-L/2 iz' 



, 2auaz{z') 



kaz{z) (24) 



L/2 {z' ~ zY + {a/2y 
that has the kernel less localized than that of Eq. ^1 



IV. COLD DEFLAGRATION EQUATIONS 

The phenomenon of cold defiagration is described by a 
collection of equations ()11|) for every magnetic molecule 
in the crystal, with the dipolar field controlling transi- 
tions being determined by the instantaneous non-uniform 
magnetization. As the full three-dimensional problem 
with a long-range interaction requires too much computer 
power, here the one-dimensional approximation will be 
made. The magnetization is considered as a function of 
the coordinate z only, i.e., the deflagration fronts are flat, 
and the dipolar field is taken along the symmetry axis as 
in Eq. (|2ip that will be used below. Of course, the dipo- 
lar field away from the symmetry axis is different, that 
will result in non-flat fronts. However, to avoid compli- 
cations in demonstrating the basic phenomenon, these 
effects will be ignored here. 

It is convenient to introduce the dimensionless time t 
and coordinate z as 



res'-; 



z/R, 



(25) 



where Frcs is the resonance relaxation rate following from 
Eq. dUl) with {S - m') / {2S) neglected. 



r — 

-L rns — 



A2 



J- m.' 



Then Eq. pT|) becomes 
d 



dt 



{z,t) ^ -F{z,t)nil,t), 



(26) 



(27) 



where F contains integral dependence on n(z, t) via 

Dzz{z,t), 



F{z,t) 



1 



En = 



2Ed 



_ (28) 
The dimensionless bias W — W/ {2Ed) , with W defined 
by Eq. ([5]), has the form 

W(z,t) = W?ext + W^''^ (z, t) = W,^, + Dzz [z. t), (29) 
where 



Ed 



{Bz — Bk 



(30) 



and Dzz{z, t) defined by Eq. (pij) can be rewritten in the 
form 



L/2 

-L/2 



l~ -\ I 1—1 2tiv(j z{z' ^t) /— ~\ 

Dzz{z,t) ^ I ^ dz J— - kaz[z,t) , 

[(P -zf + l] 



(31) 

where (7z{z, t) = 1 — 2rt_5(z, t) and L = L/ R. 

For a long sample, i — > oo, one can seek for a solution 
of Eq. (f27|) in the form of a moving front depending on 
the combined argument ^ = z — vt, where v is the front 
speed. In reduced units one has ^ = z — v*t, where the 
relation between the real and reduced front speeds has 
the form 

V = v*T,,sR. (32) 
Eq. ((TT]) for the front becomes 

v*^^FiOn, (33) 
where F(|) = (l + 4:EliW^(X)) contains 



DzziO 



d^ 



2TlV(7z{i ) 



3/2 



k<y.{0- (34) 



Eq. (|32|) makes the dependence of u on i? and Fios ob- 
vious. However, there are nontrivial parameters Ejj and 
Wcxt that enter v* . 
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FIG. 1: Propagating front of cold deflagration for Vl4xt = 2 
and Ed ~ 20. Cold deflagration starts after some ignition 
time that depends on the initial condition. 




z/(2R) 



FIG. 2: Spatial profiles of the metastable population n and 
the reduced bias W in the front for H4xt = 2 and Ed = 20. 
Everywhere in the front the system is near the resonance, 
W ~ 0. At this value of the bias periodic structures behind 
the front begin to emerge. 



For a ribbon one can introduce z = 2z/a and replace 
Eqs. ((3T|) and ((34)) by the corresponding expressions fol- 
lowing from Eq. ([M)) . 

The cold deflagration equations can be solved by dis- 
cretization that reduces them to a system of nonlinear 
ordinary first-order differential equations. The front of 
tunneling for Ejj = 20 and Wext = 2 is shown in Figs. 
[T] and [2] while more details will be given in Sec. IVII and 

lYni 



V. THE LIMIT OF STRONG DIPOLAR FIELD 

Unless a strong transverse field is applied, the numeri- 
cal value of Ed of Eq. ([28]) for Mni2 Ac is large. In par- 
ticular, in the overdamped limit Ed = Ed / {KTm') ~ lO'^ 
for to' = S" — 1. In this case one could think that F is 
negligibly small everywhere except for a very close vicin- 
ity of the resonance, so that the total relaxation and thus 



the speed of the front are very small. However, as the 
numerical solution shows, the system finds the way to 
relax faster by forming a front region of a width I ^ R 
where W is small and F is of order one. In this extended 
region resonant tunneling transitions take place. Beyond 
the front region W deviates from zero and F becomes 
negligibly small. As a result, n{^) changes practically 
only within the front core. 

Basing on these insights, one can construct a pertur- 
bative expansion in powers of 1/Ed and show that the 
solution n(^) and the front speed v* become independent 
of Ed for Ed ^ 1. One can search for az{£,) in the form 



(0 



r(l) 



(35) 



and similarly for n(^) 



1 - <^z{0] /2. The term ^^"^(C) 
is defined by the condition that within the front region 
—1<£,<1, with the width I to be determined self- 
consistently, the contribution of a^f^ (^) into the bias W 
IS zero. If this is fulfilled, the term AEj^W'^ in the de- 
nominator of Eq. ([^5]) is of order one due to the cor- 
rection a')}'\C)-, so that in the front region F is of or- 
der one. In the region before the front, I < ^, one has 
~ 1 and (J^^\s,) = Uzi = — 1. Everywhere behind 
the front, ^ < — one has final values n^°)(|) = Uf and 
of\C) = cFzf that are to be determined. In the front 
region the condition W — with Eq. yields the in- 
tegral equation 



(36) 



where 



H^(0)(^) = We: 



~l ~N 2 

s, -{] +1 



-1 3/2 



(37) 



This equation determines the zero-order profile ai^\^), 
including I and CTj/. 

Eq. ([36]) can be solved numerically by discretizing the 
integral using -I- 1 equidistant points within the inter- 
val (-^,^ given by I, = - T+ 2li/N, i = 0,1,..., N. 
The value at the right end of the interval is fixed by the 
boundary condition af \£,j^) = (J^i = — 1. Thus there are 
total A^ + 1 unknowns including I, that can be found by 
solving the system of A^ -I- 1 equations VF^"' (^J = with 
i — 0,1, . . . , N. Note that this system of equations is non- 
linear because of L In this way one finds the zero-order 
magnetization profile in the front and the magnetization 
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behind the front a^f for any Woxt > 0. In particular, for 
Mni2 Ac one obtains I = 0.848. 

On the other hand, one can find important analytical 
results if one searches for the solution in the form 



-/(C), 



(38) 



where f{il) — ±1. Substituting this into Eq. (|37p one 
obtains the equation for /(^) 



e +1 



-1 3/2 



-TTV {a,f 



where oi'^'^'' = ^t^v — k. One can see that there are even 
and odd terms in ^ in this equation and f{$) is odd. 
The even and odd parts of this equation should turn to 
zero independently of each other. For the even part one 
obtains the equation 




(39) 



= We. 



zz 2 



(40) 



that with a , 



1 yields CT./ = 1 - 2Wcxt/^i7'^ and 

1 - Ozj _ VFext 



D 



(cyi) 



(41) 



for the fraction of metastable molecules behind the front. 
Note that since < n/ < 1, this solution only exists for 

0< Wext <^i7'^. (42) 
The odd part of Eq. (IHS)) yields the equation 



(r-i)= 



1 



3/2 27ri/ 



12 



(43) 



that defines /(C) and I. They can be found numerically 
by discretization as described above. The expression for 
n'^°^(^) in terms of /(^) following from Eq. ([55)1 has the 
form 



r.(0)(e) 



1 



1 



1 



2 \ ' f)('=yi) 



1 



2\' Dit\ 



no- (44) 



Using the method of Ref. 0, one can show that the ap- 



proximate solution for /(^) valid for 

/app(0 - ill app, 



<C 1 has the form 
(45) 
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-0.5 



Ml 





y thick slab 
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FIG. 3: Normalized magnetization profile in the cold deflagra- 
tion front for long Mni2 Ac crystals of cylindrical and thick- 
slab shape. Approximate result for the cylinder is shown by 
a dashed line 
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FIG. 4: Normalized energy bias W in the cold deflagration 
front in long Mni2 Ac crystals in the limit of strong dipolar 
fleld. In the central part of the front VF = so that resonance 
transitions take place. 



where 



^app — 



(46) 



For a Mni2 Ac cylinder one obtains Zapp = 0.7137. 

Numerically found /(^) and its approximation /app(C) 
for a Mni2 Ac cylinder is shown in Fig. [3] together with 
/(^) for a thick slab (ribbon) discussed at the end of this 
section. 

The speed of the front v* can be found by considering 
the effect of the correction o'i^''(^), although one does 
not need to evaluate this correction explicitly. In Eq. 
([55)1 one has to keep cri^^(^) in F, because otherwise the 
denominator turns to zero. This makes F zero order in 
\/Ed- On the other hand, dn/d^ and n outside F can be 
taken at zero order in l/Eo- One thus can rewrite Eq. 
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(155)1 in the interval — Z < ^ < Z in the form 



9(0 



(47) 



There is a point inside the interval —l<^ < I where W 
changes its sign. At this point 1/F reaches its minimal 
value 1. On the other hand, this point can be deter- 
mined as the minimum of the rhs of this equation. Then, 
obviously, 



5(0 



Using Eq. (j44|) one obtains 
One can see that for Wcxt - 



Q 



D 



(cyl) 



one has min 



(48) 



(49) 



0, 



achieved sA — I where / = — 1. This yields v* = at 
W^ext = 0. In the general case one has to investigate 



5(0 



1 



Q + 1 + fiO 



(50) 



no 

Since g'{-T) = at 

Q = Qc - r (-0//"(-0, (51) 

one concludes that for Q < Qc the minimum is achieved 
at ^ = — Z and thus 



n-i) 



W^cxt f'i-l) 



(52) 



For a cylinder of Mni2 Ac one has 2/ f'{—l) — 2.31 and 
Qc — 0.809. According to Eq. (^51) . the latter translates 
into W^ext,c = 3.03. Then Eq. jSHl) yields the value of 
the corresponding bias field B^^c — Bzk = 15 mT. For 
Bzk ^ Bz < Bz^c, the front speed in real units obtained 
with the help of Eqs. and (HOI) is given by 

Bz — Bzk 2 



V = RTrc 



ijz — Dz 



(53) 



Bzk /'(-O 

For Qc < Q (and thus Woxt,c < Wcxt) one has to find 
min g{^) from the condition g'{S,) — that leads to 
somewhat smaller front speeds than given by the for- 
mulas above. However, the laminar regime of the cold 
deflagration fronts breaks down at the external bias 
smaller than Woxt,c, so that the results of this section 

for Wext,c < Woxt are irrelevant. 

Let now consider the slab geometry. From Eq. ([M]) 
with z = 2z/a instead of Eq. ([57]) one obtains the equa- 
tion 



— h Arctan ( ^ 



(1 



Arctan i^ + l 



(54) 
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Discretization: 
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Discretization: 
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o A'=400 



FIG. 5: Reduced front speed v* vs the reduced bias VI4xt for 
different discretizations. In (a) for Ed = 100 and Wcxt < 1 
the numerical results are in a good accordance with the large- 
Ed formula (straight line). 



The numerically obtained result for the normalized mag- 
netization profile /(^) is shown of Fig. [31 compared to 
that of a cylinder. Since the kernel in the integral equa- 
tion for the slab is less localized for a thick slab than for 
a cylinder, the front width I — 1.61 for a thick slab is 
larger than I = 0.848 for the cylinder. All formulas ob- 
tained above are valid for a thick slab as well, however 
with different constants: 2/f'{-J) = 4.79, Qc = 0.540, 
and W?ext,c = 2.24. 



VI. NUMERICAL RESULTS 



As mentioned at the end of Sec. IIVI the cold- 
deflagration equations can be solved by discretization re- 
ducing them to a system of ordinary differential equation. 
Numerical calculations use the semi-infinite geometry in- 
cluding the region of length —L/2 < z < L/2 where 
equations are solved plus the range L/2 < z < oo where 
the magnetization is fixed to a z — corresponding to 
the metastable state. The latter is needed to create the 
dipolar field in the main region — _L/2 < z < L/2 that 
corresponds to the semi-infinite sample. This allows to 
operate on shorter samples that saves computer time. 
When the deflagration front reaches z — L/2, it cannot 
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go further, so the results near this point become unphys- 
ical and should be ignored. 

The first thing revealed by computations is that for 
large values of Eo it is very important to prepare the 
system in the initial state close to the actual front, with 
W ^ within the front core. The further is the initial 
state from this optimal state, the more time (ignition 
time) the system needs to adjust so the the front could 
start moving across the sample. For initial states far 
from the front states, the ignition time can be so long 
that there is a significant off-resonance relaxation in the 
bulk of the crystal during it. For smaller dipolar fields 
such as Eu ~ 3 (that can be achieved by applying a 
strong transverse field to increase A) ignition of the fronts 
is much easier. Another possibility to ignite the cold 
deflagration is to slowly sweep the external field in the 
positive direction, approaching the resonance^i^ that will 
be considered later on. 

Computations for large Eo and not too strong bias 
Woxt corroborate semi-analytical results of the preceding 
section. For z not too close to the ends of the interval 
—L/2 < z < L/2, the variables indeed depend on ^ = 
z ~ vt, as it should be in a moving front. 3D plots of 
n{z,t) are smooth and look qualitatively as in Fig. [TJ 
and the ignition time can be reduced to zero by a better 
choice if the initial state. The metastable population 
n{t) averaged over the length of the crystal is almost flat 
during the ignition time, then it decreases linearly as the 
front travels through the crystal, then becomes nearly 
flat again after the front arrives at the right end of the 
interval -L/2 < z < L/2, see Fig. 16 of Ref. [U for the 
standard magnetic deflagration and Fig. 4 of Ref. [l^ for 
the cold deflagration. The front speed can be obtained 

^ — L I (^arrival ^ignition) ■ 

The reduced front speed i;*, Eq. ([5^ . vs Woxt, Eq. 
is shown for Eq = 100 in Fig. Eli. For not too 
large bias, Wcxt ^ 1, the numerical results are in a good 
agreement with Eq. (|5^ shown by a solid line. In this 
region the numerical results do not depend on the number 
of discrete points used in the solution. Similar results for 
Ee, = 20 in Fig. [5)3 are further from the theoretical curve 
(not shown) because the condition of a strong dipolar 
field is not fully satisfied. Also there is some nonzero 
speed in the region Wcxt < that is, however, quickly 
decreasing with the negative bias. 

With increasing the external bias, the laminar solution 
in the form of a smooth moving front loses its stability. 
In Fig. [2] one can already see wiggles behind the front 
that represent frozen-in spatial structures with the period 
of order R. With increasing Woxt or Ed, these features 
progress and the region of transitions moves with oscillat- 
ing speed, see Fig.[Bl To distinguish this transition region 
from the true front, it was called "wall of transitions":^^ 
It would cost significant additional efforts to find out an- 
alytically or numerically whether the transition from the 
laminar to non-laminar regime with increasing Woxt is 
gradual or there is a threshold. One important obser- 
vation is that the spatial structures behind the front are 




£n= 100 10000 



FIG. 6: Propagating front of cold deflagration for Ed = 100 
and Wext = 5. Here ignition time was eliminated by good 
choice of the initial condition. The front speed is oscillating 
and there are spatially-periodic structures behind the front. 
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FIG. 7: Front speed v* vs Ed in the non-laminar regime for 
Wcxt = 5. 



discontinuous on z, while the analytical result of Eq. ([48]) 
was obtained based on the assumption that the solution 
is continuous. 

As the laminar regime breaks down, the instability is 
manifested by the dependence of the results for v* in Fig. 
[5| on the number of discretization points. With increas- 
ing A'' the discontinuities in w*(Woxt) become smaller but, 
unfortunately, increasing N is limited by computing re- 
sourses. For larger Woxt the front speed reaches a plateau 
that depends on Eo , and curves with different discretiza- 
tions converge again. 

In the non-laminar regime the magnetization in the 
front cannot completely adjust so that bias in the front 
core would be very close to zero and the resonance transi- 
tions could occur at a rate close to the maximal rate Fj-es- 
This is the reason why v*{Wext) drops after reaching a 
maximum, as the instability begins. Still the very exis- 
tence of the front in this case suggests that the system is 
closer to the resonance in this region than in the others. 
The values of W in the front should be of order 1, so that 
in Eq. ([28]) one has F ~ ^/E-^- This is supported by the 
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FIG. 8: Wall of tunneling at Ed = 100, ignited by slow sweep 
of the bias field, e = 30. The process starts at M4xt ~ 5. 
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FIG. 9: Dipolar field in the crystal at different values of Woxt 
during sweep. After Wext = 5 the front starts moving and 
dipolar field in the sample becomes discontinuous. 

computations shown in Fig. [7| In the plateau region (in 
particular for Wext — 5) the front speed fits to 

V* ~ 8/i?|) (55) 

for large Eu- 

VII. COLD DEFLAGRATION INITIATED BY 
FIELD SWEEP 

As was mentioned above, for most of initial conditions 
the development of the cold deflagration front requires a 
very long ignition time. If the initial condition is chosen 
in a special way close to the actual front, the process 
starts immediately. However, one cannot find a practical 
way to prepare such initial state. 

Fortunately, as was found in Ref. [l^ the front can be 
ignited by a slow time-linear sweep VFcxt — vwt start- 
ing with a value of Wcxt that ensures W <Q everywhere 
in the sample, see Fig. [5] The sweep rate can be con- 
veniently parametrized by e = ttA^/ {2hvw) and slow 
sweep requires e ^ 1. As Wcxt increases, the condition 



W = Q would be first reached at the end of the sam- 
ple, then the resonance point would move into its depth. 
However, transitions induced by the sweep (seen in Fig. 
|S] before the ignition) change the dipolar field so that 
the system does not cross the resonance in the region 
near the end of the sample, although it becomes close to 
the resonance in the increasingly broad region, see Fig. (S] 
The reason for this is that flipping spins in a small region 
near the end of the sample do not significantly change the 
dipolar field from the surface of the crystal, the integral 
term of Eq. (|2T|) . but strongly change the local contri- 
bution, last term in this formula. Increasing az due to 
resonance transitions leads to the decrease of the local 
term that creates a negative dipolar bias and prevents 
the system from crossing the resonance. After some time 
the region close to the resonance becomes broad that is 
similar to the structure of the cold deflagration front, see 
Fig. [51 In this way the initial state for the cold defla- 
gration is being prepared. The front starts as the bias 
reaches the "magic" value of VFoxt that weakly depends 
on Ed- For En = 20 one has Wext — 4.3 (that corre- 
sponds to Bz — Bzk — 19 mT) and for Ejj — 100 one has 
Wcxt — 5. At such a strong bias there is a quasiperiodic 
spatial structure with discontinuous magnetization and 
the dipolar field behind the front, as shown in Figs. [5] 
and[5] One can see in Fig.[S]that in the moving front the 
bias is slightly below zero. This means that the system is 
somewhat off-resonance and this is the reason for a small 
front speed in this regime, as shown in Fig. [5] 

The next question is how to ignite cold deflagration for 
arbitrary values of the external bias Wcxt given by Eq. 
(150)) . The answer is to sweep up to this value of Wcxt 
and then to stop this (global) sweep at some imax.o- After 
that further sweep Bz locally near the end of the sample, 
z — ~L/2, using a small coil. For the coil of radius Rq 
placed at zq (the axis of the coil coincides with the axis 
of the cylinder) the local addition to Wext can be written 
in the form 

<5Wcxt(2, t) = — -2 -^Vw [t - imax,o) (56) 

i?^ + (z - zo) 

with Z{) w — L/2. Numerical calculations with Rq — R 
and zq = —L/2 show that, indeed, with this method one 
can ignite fronts at different biases Wext, including the 
laminar regime, see Fig. llOl Here the front is much faster 
than in Fig. [51 in accordance with the results for the front 
speed in Fig. [5^. Total bias in the sample at different 
times is shown in Fig. [TlJ For instance, r = Fics^ = 5000 
corresponds to the stage of the global sweep. All other 
times correspond to the local-sweep stage, since the bias 
curves converge on the right side of the sample where 
<^Wcxt is small. Local sweep near the left end creates 
an initial state for the front to start, as the bias curves 
are approaching zero in a progressively large region (blue 
curves). As the front starts moving (red curves), the bias 
becomes positive on the left with non-laminar features 
near the end. But in the depth of the sample the front 
is laminar corresponding to Wcxt = 1- 
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FIG. 10: Laminar front of tunneling at En = 100 and M4xt = 
1 ignited by slow local sweep of the bias field, (a) Overview; 
(b) Zoom of the front region. 
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FIG. 11: Total_bias at differeril^ times in the laminar front of 
tunneling at Ed ~ 100 and Wcxt ~ 1 ignited by slow local 
sweep of the bias field. Blue lines correspond to times of 
adjustment before ignition and red lines correspond to the 
moving front. 



VIII. DISCUSSION 

In the main part of the paper it was shown that elon- 
gated crystals of molecular magnets (practically Mni2 
Ac) can exhibit moving fronts of dissipative spin tunnel- 
ing at biased resonances, Eq. ^ with fc > 0. This phe- 



nomenon is resembling magnetic deflagration , ^"i^-'^ only 
the relaxation rate is controlled not by the temperature 
but by the dipolar field evolving self-consistently and 
bringing the spins in the front region on and off tun- 
neling resonance. Like deflagration, it leads to destruc- 
tion of the initial metastable ordered state (however in 
general incomplete), this is why it can be called "cold 
deflagration" . 

Of course, transitions at biased resonances result in en- 
ergy release and warming of the sample, so that the two 
mechanisms can coexist. In fact, magnetic deflagration 
was observed in crystals of Mni2 Ac thermally isolated 
so that the warming of the sample is efficient. With- 
out thermal isolation, cold deflagration does not face this 
competition. To further reduce heating, it is preferable 
to work at low bias, such as fc = 1. 

There are two regimes of cold deflagration: laminar 
regime at low bias — Bzk and non-laminar regime at 
high bias. In the laminar regime the magnetization in the 
front adjusts so that the dipolar fleld B'^^^ together with 
the external fleld B-^ creates a nearly zero bias for the 
resonance transitions in the front region with the width 
of order R, the transverse size of the crystal. In the 
laminar regime the magnetization and dipolar fleld in 
the sample are continuous and both the front speed and 
the magnetization (metastable population) behind the 
front can be found analytically in the practical limit of 
the strong dipolar field, Eqs. (|53)) and (|4T|) . Remarcably, 
both of these quantities do not depend on the strength 
of the dipolar field En in this region. 

In the laminar regime the estimation for the front 
speed is u ~ RTres, where T^es is the transition rate 
at resonance, W = Q in Eq. p^ . At the boundary be- 
tween the over- and underdamped regimes, A '-^ T^' and 
thus Tres ~ Tm' (that is realized in the transverse field 
3 T in Mni2 Ac at the k = 1 resonance) cold defiagra- 
tion already beats the regular "hot" defiagration. Indeed, 
the latter has the speed v ~ lT{Tf ), where I depends on 
the thermal diffusivity but experimentally is comparable 
with R and r(T/) is the thermal activation rate over the 
barrier at the fiame temperature. Since r(T) at high 
temperatures is determined by the rates of transitions 
between adjacent levels near the top of the barrier that 
are smaller than the same for low-lying levels such as r„j' 
(one has r(T/) - 10^ s^^ and Ts-i ^ l(f s^^) the hot 
deflagration loses. 

At higher bias B^ — B^k the laminar regime breaks 
down, the dipolar fleld cannot fully adjust to provide a 
nearly zero bias in the front's core, and the magnetiza- 
tion and the dipolar fleld become discontinuous. There 
are frozen-in quasiperiodic spatial structures behind the 
front. Accordingly, the front speed dramatically drops, 
see Fig. [51 especially in the case of a strong dipolar field. 
There is no analytical solution in this range but the fit to 
the numerical results yields v cc The boundary be- 

tween the laminar and non-laminar regimes corresponds 
to Bz - Bzk = 5-10 mT. 

It was shown that cold defiagration can be ignited by 
the local sweep of the field B^ near an end of the crystal. 
This local field can be produced by a small coil with 
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increasing current. 

Another condition for the observability of the cold de- 
flagration is sufficiently strong transverse field that allows 
tunneling transitions via low-lying levels. 

At nonzero temperatures the rate of cold deflagration 
should increase because of the activation to higher levels 
providing a higher transition probability, see Eq. (10) of 
Ref. M 



Acknowledgements 



The author is indebted to E. M. Chudnovsky for many 
stimulating discussions. 

This work has been supported by the NSF Grant No. 
DMR-0703639. 



1 T. Lis, Acta Crystallogr. B 36, 2042 (1980). 

^ A. C. R. Sessoli, D. Gatteschi and M. A. Novak, Nature 

(London) 365, 141 (1993). 
^ J. R. Friedman, M. P. Sarachik, J. Tejada, and R. Ziolo, 

Phys. Rev. Lett. 76, 3830 (1996). 

* J. M. Hernandez, X. X. Zhang, F. Luis, J. Bartolome, 
J. Tejada, and R. Ziolo, Europhys. Lett. 35, 301 (1996). 

^ L. Thomas, F. Lionti, R. Ballou, D. Gatteschi, R. Sessoli, 

and B. Barbara, Nature 383, 145 (1996). 
® A. Morello, E. L. Mettes, F. Luis, J. F. Fernandez, 

J. Krzystek, G. Aromi, G. Christou, and L. J. de Jongh, 

Phys. Rev. Lett. 90, 017206 (2003). 
^ D. A. Garanin and E. M. Chudnovsky, Phys. Rev. B 78, 

174425 (2008). 

* S. McHugh, R. Jaafar, M. P. Sarachik, Y. Myasoedov, H. 
Shtrikman, E. Zeldov, R. Bagai, and G. Christou, Phys. 
Rev. B 79, 052404 (2009). 

^ N. V. Prokof'ev and P. C. E. Stamp, Phys. Rev. Lett. 80, 
5794 (1998). 

^° A. Cuccoli, A. Fort, A. Rettori, E. Adam, and J. Villain, 
Eur. Phys. J. B 12, 39 (1999). 

J. J. Alonso and J. F. Fernandez, Phys. Rev. Lett. 87, 
097205 (2001). 

J. F. Fernandez and J. J. Alonso, Phys. Rev. Lett. 91, 
047202 (2003). 

P. C. E. Stamp and L S. Tupitsyn, Phys. Rev. B 69, 014401 
(2004). 



L S. Tupitsyn, P. C. E. Stamp, and N. V. Prokof'ev, Phys. 
Rev. B 69, 132406 (2004). 

J. F. Fernandez and J. J. Alonso, Phys. Rev. B 69, 024411 

(2004) . 

J. F. Fernandez and J. J. Alonso, Phys. Rev. B 72, 094431 

(2005) . 

W. Wernsdorfer, T. Ohm, C. Sangregorio, R. Sessoli, D. 
Mailly, and C. Paulsen, Phys. Rev. Lett. 82, 3903 (1999). 

D. A. Garanin and E. M. Chudnovsky, Phys. Rev. Lett. 
78, 097206 (2009). 

Nurit Avraham, Ady Stern, Yoko Suzuki, K. M. Mertes, 
M. P. Sarachik, E. Zeldov, Y. Myasoedov, H. Shtrikman,l 

E. M. Rumberger, D. N. Hendrickson, N. E. Chakov, and 
G. Christou, Phys. Rev. B 72, 144428 (2005). 

Y. Suzuki, M. P. Sarachik, E. M. Chudnovsky, S. 
McHugh, R. Gonzalez-Rubio, N. Avraham, Y. Myasoedov, 
E. Zeldov, H. Shtrikman, N. E. Chakov and G. Christou, 
Phys. Rev. Lett. 95, 147201 (2005). 

D. A. Garanin and E. M. Chudnovsky, Phys. Rev. B 76, 
054410 (2007). 

D. A. Garanin and E. M. Chudnovsky, Phys. Rev. B 56, 
11102 (1997). 

D. A. Garanin, arXiv;0805.0391 (2008). 

E. M. Chudnovsky, D. A. Garanin, and R. Schilling, Phys. 
Rev. B 72, 094426 (2005). 



